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Abstract 

A numerically exact path integral treatment of the absorption and emission spectra of open quantum 
systems is presented that requires only the straightforward solution of a stochastic differential equation. 
The approach converges rapidly enabling the calculation of spectra of large excitonic systems across the 
complete range of system parameters and for arbitrary bath spectral densities. With the numerically exact 
absorption and emission operators one can also immediately compute energy transfer rates using the multi¬ 
chromophoric Forster resonant energy transfer formalism. Benchmark calculations on the emission spectra 
of two level systems are presented demonstrating the efficacy of the stochastic approach. This is followed 
by calculations of the energy transfer rates between two weakly coupled dimer systems as a function of 
temperature and system-bath coupling strength. It is shown that the recently developed hybrid cumulant 
expansion is the only perturbative method capable of generating uniformly reliable energy transfer rates 
and spectra across a broad range of system parameters. 
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I. INTRODUCTION 


The far-held absorption and emission spectra are standard experimental tools in the character¬ 
ization of excitonic systems. The temperature and solvent dependence of these spectra are often 
used to extract a wealth of information on, for example, the microscopic geometry of the con¬ 
stituent chromophores, the coupling strength between the excitonic system and its environment, 
as well as the relative importance of heterogeneous broadening mechanisms.—”— Despite this wide 
applicability, at present there are still relatively few theoretical approaches that are capable of pro¬ 
viding uniformly reliable estimates of the spectra of open quantum systems. Only in the limiting 
case that the spectrum arises from a single isolated electronic transition can the exact absorption 
and emission spectrum be obtained analytically (up to a numerical integration) through cumulant 
expansion techniques^ However, in the more common setting, wherein the excitonic complexes 
are comprised of multiple coupled chromophores, then one must, in general, resort to numerical 
methods. Unfortunately, there is no numerically exact approach currently available for systems 
containing more than a few chromophores that is valid over a large range of the parameter space. 
As a result, comparisons between many interesting experimental systems and their corresponding 
microscopic theoretical models are often out of reach. One of the central aims of this work is to 
fill this gap. Here we present an efficient path integral approach that allows one to compute the 
numerically exact absorption and emission spectra of multi-chromophoric open quantum systems. 

Due to the lack of robust exact methods, one often turns to perturbative techniques. As 
detailed in the preceding papers of this series (henceforth referred to as papers and II^) many 
of the standard approximate approaches are capable of generating reliable absorption spectra, as 
only the short time dynamics of a factorized system-environment initial state are required. The 
emission spectrum, however, presents a much more challenging problem. In this case, the real¬ 
time dynamics evolve from the correlated equilibrium state of the entire excitonic system and its 
environment. Unless the system-bath coupling is very weak, perturbative treatments often generate 
qualitatively incorrect emission spectra, and generally become even worse as the temperature is 
loweredj^i^ One of the major results of the previous papers in this series was a hybrid perturbative 
method capable of providing reliable emission spectrum over a broad range of system parametersi^ 
The approach is not fully perturbative in that it combines the knowledge of the numerically exact 
equilibrium reduced density matrix -which can be obtained relatively easily through imaginary 
time path integral methods-^ with an approximate cumulant expansion of the remaining real time 
dynamics. This initial state correction becomes essential at low temperatures or strong coupling. 
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The hybrid cumulant expansion (HCE) thus greatly extends the parameter regimes accessible to 
perturbative methods and generally improves the quality of the results. 

In the context of numerically exact treatments of the emission spectrum, such as the hierarchy 
equation of motion (HEOM) or quasi-adiabatic path integral approach (QUAPI),— the problem of 
a correlated system-bath initial state is overcome by simply preparing a factorized state sufficiently 
far in the past such that the system has reached equilibrium at time zero. As a result, these 
approaches require an initial lengthy propagation of the reduced density matrix to equilibrium 
before the dynamics of the dipole correlation function can be calculated.-! Furthermore, as these 
approaches rely on efficient representations of the inffuence functional, they are generally restricted 
to environments that are not strongly non-Markovian. The main result of this paper is a stochastic 
path integral approach that circumvents many of the restrictions imposed by other numerically 
exact methods, and in particular, is applicable for arbitrary spectral densities and temperatures. 

While an exact calculation of the absorption and emission spectra is important in its own right, 
it also provides an additional benefit. That is, one can immediately compute energy transfer 
rates between weakly coupled excitonic systems using the multi-chromophoric Forster resonant 
energy transfer (MCFT) formalism. The MCFT framework is a generalization of the standard 
Forster theory to the situation where the donor or acceptor complex consists of multiple coupled 
chromophores,— and has gained recent attention as this scenario appears to be one of most com¬ 
mon motifs employed in the highly efficient energy transfer networks found in biological systems. 
For example, the light harvesting systems found in both green and purple bacteria are comprised of 
independent complexes of strongly coupled chromophores that form the base units for large-scale 
energy transfer networks.— 

In Sec. m the details of the MCFT formalism are presented. There it becomes apparent that 
the key quantities necessary for computing energy transfer rates are generalized operators related 
to the absorption spectrum of the acceptor complex and the emission spectrum of the donor. Then 
we present the path integral treatment of the absorption operator and demonstrate that it may 
be efficiently computed as the solution to a straightforward stochastic differential equation. This 
approach is then generalized to the emission spectrum by taking advantage of the detailed balance 
condition that relates the emission operator to its corresponding absorption operator evolving in a 
complex time. Following these formal developments, numerical calculations are presented for model 
two level systems that can be reliably benchmarked against the HEOM approach. The temperature 
dependence of the emission spectrum is presented, followed by systematic calculations of the MCFT 
rate as a function of the temperature and system-bath coupling strength. It is observed that the 
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hybrid cumulant expansion techniqne developed in paper II is the only perturbative approach 
that provides uniformly reliable results for the energy transfer rates^^ In a forthcoming work, the 
path integral and HCE methods will be used to provide the first systematic analysis of the energy 
transfer rates between two B850 complexes in the light harvesting system LH2.— 

II. MCFT FORMALISM 

The MCFT formalism has been expounded in the previous papers in this series. Here we 
provide only the salient details necessary to keep the presentation self contained. The total system 
is composed of a donor complex consisting of chromophores that is weakly coupled to an 
acceptor complex of chromophores. The Hamiltonian for the entire donor-acceptor system is 
then 


H = H^ + H^ + , ( 1 ) 

where denotes the Hamiltonian operator of the multi-chromophoric donor (acceptor) com¬ 

plex along with its associated thermal environment. The excitonic coupling between the donor and 
acceptor systems is characterized by which, within the local basis of the single-excitation 

subspace of the donor and acceptor, is given by 

No Na 

n=l m=l 

The Hamiltonian of an individual complex is modeled as a general open quantum system, 

H- = H^ + H^ + H2, (3) 

where the label a E {D, A) serves to distinguish between the donor and acceptor systems. The free 
excitonic Hamiltonian of each complex is given by 

iVc ATc 

Hf, = ^ ^ (e^ -|- \arn) (o^ml + ^ ^ ^nm Ic^n) {otm\ i (4) 

m=l n^m 

where Em is the excitation energy of the m-th chromophore, trim denotes the intra-complex elec¬ 
tronic couplings and Am is the environment-induced reorganization energy. The free bath Hamil¬ 
tonian is 

No, 

= ( 5 ) 

m=l k 
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where denotes the respective creation (annihilation) operator of the k-th mode of the 

bath with frequency 10 “^ and coupled to chromophore m on the excitonic complex labeled by a. 
The system bath coupling is linear in the bath coordinates, and assumed to modulate only the 
excitation energies, 

HI = E E [btk + b^,k) , (6) 

m=l k 

where = \am} (oml and denotes the coupling strength. 

Assuming that the exciton lifetime is much longer than the timescale associated with the energy 
transfer, then relaxation to the ground state can be safely ignored, and the population transfer 
rate between the donor and acceptor systems is given by the MCFT rate formula, 


k = 2Re 


f 


dt Tr 


{H 




(7) 


which can be easily obtained from the golden rule expression as shown in section II of paper L— The 
absorption operator of the acceptor, and emission operator of the donor, E^(t), appearing 

in Eq. [7]are formally defined as: 


l^{t) = TVb 
E°(t) = Trb 




( 8 ) 

(9) 


In the case of the absorption operator, the initial density matrix corresponds to a factorized state of 
the system and bath, /?^ = /g 0 l^hi assumption of a Franck-Condon transition 

from the ground state. The steady-state emission, however, occurs after the total system has 
equilibrated within the single excitation manifold. Thus, the initial state in Eq. 0 corresponds to 


IS 


the equilibrium state of the entire system and bath, where = Tr 

the partition function of the donor. As discussed extensively in the previous papers of the series, 
the difference in the initial states is the key feature that distinguishes the absorption from the 
emission operator, with the correlated initial condition in Eq. [9] leading to a substantially more 
involved calculation. 


A. Detailed Balance 


The absorption and emission spectra obey a well-known detailed balance condition, and it is 
readily apparent that their corresponding operators in Eqns. [8]and[9]must obey a similar relation. 
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In the frequency domain, the detailed balance condition for the operators reads 


. ( 10 ) 

where Z = jZ^. Thus, in principle, knowledge of the absorption operator allows for a straight¬ 
forward determination of the corresponding emission operator. However, in practice, the thermal 
prefactor exponentially amplifies any error in the absorption data leading to an ill-conditioned 
numerical problem. As a result, Eq. (TOlis generally of little practical use outside of the very high 
temperature limit. 

An alternative approach can be based on the observation that in the time domain, the detailed- 
balance condition takes the form 

= ( 11 ) 

where the asterisk denotes complex conjugation. That is, through the straightforward substitution, 
t —?• t — ih/3, the time evolution of the emission operator of the donor is equivalent to that of the 
absorption operator, except that the dynamics evolves in complex time rather than in purely real 
time. In contrast to the frequency-domain detailed balance relation, the time-domain version in 
Eq. dUis free from numerical instabilities and forms the basis for the developments presented here. 
Here we employ the path integral formalism to develop an exact and efficient numerical treatment of 
the spectral operators rather than pursue perturbative approaches as were explored in the previous 
papers of this series. In the following subsection, the stochastic path integral representation for the 
absorption operator is presented, and then generalized to the case of emission through the rotation 
from real time to complex time suggested by Eq. fTTl 

B. Absorption Operator 

As can be seen from Eq. [HI the absorption operator does not require the full time evolution 
of the reduced density matrix. The bath evolves both forward and backward in time, but the 
system is only propagated forward in time. As a result, we can still take advantage of the influence 
functional formalism from the path integral approach to open quantum systems,— but only require 
a single path for the system variables. Thus the absorption operator can be determined from the 
path integral expression 

U^{t) = j V[a]e^^^^^^F[a] , (12) 
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where 5',^[o'] denotes the action associated with the free system Hamiltonian of the acceptor, 
and the standard boundary conditions of the paths have been suppressed for clarity. The 
Feynmann-Vernon influence functional, T[(t], obtained by integrating out each of the Va indepen¬ 
dent baths is given by 


Na 

F[o-] = exp 

n=l 






(13) 


All of the microscopic details of the baths that are relevant to the system dynamics enter through 
their respective correlation functions in the influence functional, which take the standard form, 


1 /■°° 

Cn{t) = — / dio Jn{io) [coth.{h/3io/2) cos{iot) — ism{iot)] , (14) 

TT Jo 

with the spectral density function, 

Jn{^) = - U)n,k) ■ (15) 

2 ^ u:n,k 
k ’ 

One of the great benefits of the path integral formalism is that it places no restrictions upon the 
functional form of the spectral density as opposed to many other approaches to open quantum 
systems. 

Following our previous developments on the equilibrium reduced density matrix,— as well as those 
of several others on the full real time dynamics of the density matrix,— the nonlocality present 
in the influence functional can be substituted for local interactions with stochastic auxiliary fields, 
which can then be efficiently sampled through Monte Carlo methods. Formally, this is affected by 
applying a separate Hubbard-Stratonovich transformation to each of the Va terms in the influence 
functional. Then Eq. [13] can be exactly rewritten as 

^H=n J Jq ^ ^ 

(16) 


where Wn represents the normalization constant of the Gaussian functional integral associated with 
the n-th bath. The path integral involving the system variables is now completely local in time and 
the auxiliary fields can be reinterpreted as a source of colored noise driving the system dynamics. 
Thus individual samples of the absorption operator can be simply and straightforwardly calculated 
through the solution of a stochastic differential equation, 

, ( 17 ) 
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subject to the initial condition /O^(0) = Ig- The stochastic Hamiltonian is given by, 


= + (18) 
n 

and the scalar, complex-valued, Gaussian noise terms obey the correlations, 

iUt)) = 0 , 

{^n{t)U{t')) = 6nmCn{t - t!)/% . (19) 

The exact time evolution of the absorption operator is obtained after averaging the stochastic 
dynamics over realizations of the noise, I"^(t) = Here we have assumed that each of 

the baths is independent of all others. To include correlated baths, one need only to replace the 
delta function in Eq. [19] with the desired spatial correlations. The generation of complex, Gaussian 
colored noise is discussed in Appendix O 

The stochastic path integral equation for the absorption operator bears some similarity to the 
non-Markovian quantum state diffusion (NMQSD) approach recently proposed to compute the 
zero temperature absorption spectrum of excitonic systems.— While NMQSD is formally exact, 
all practical implementations to date have relied on approximations. In contrast, the present 
Eq. [T3 is both formally and numerically exact, although it may prove fruitful to further explore 
the connections between the two approaches. 


C. Emission Operator 


Due to the correlated initial state, the calculation of the emission operator is considerably more 
involved than that of the absorption. The propagator in the path integral representation is,— 

U^{t,np) = j V[a\ J V[a]e-^^oW^-^^o’°l^]F[a,a] , (20) 

E D 

where Sq ’ denotes the Euclidean action of the donor system Hamiltonian associated with the 
initial imaginary time propagation to the equilibrium state. The influence functional now con¬ 
tains three contributions from the respective propagations in real time, imaginary time, and the 
correlations between the two, 



The bath correlation function is defined for complex arguments through the analytic continuation 
of Eq. [13] asr^ 



( 22 ) 


where z = t — ir, and 0 < r < fij3. 

As is readily seen the path integral expression for the emission operator is considerably more 
complicated than the corresponding result for the absorption operator. Additionally, the coupling 
between the real and imaginary time paths in the influence functional prevents a straightforward 
application of the Hubbard-Stratonovich transformation as was used previously for the absorption 
operator. Fortunately, a simplification is possible. The detailed balance relation in Eq. [TT] suggests 
that the emission operator may be computed in an identical manner to the absorption through the 
introduction of a complex time variable z = t — ifil3. Indeed, with this substitution in the path 
integral expressions above, the emission propagator may then be defined in an analogous fashion 
to Eq. [12] except along a time-ordered contour in the complex time plane,— 



(23) 


The influence functional similarly simplifies to. 



dz" ED(a(z'))ED(cT(z"))C„(z' - z") . (24) 


With these results, the Hubbard-Stratonovich transformations and stochastic Schrodinger equa¬ 


tion are then formally equivalent to those of the absorption operator presented in Sec. IIIBl The 


complex time evolution of samples of the emission operator obeys the equation 


dz n 


(25) 


where the integration proceeds from time z = 0 to the complex time z = t — ihl3 subject to the 
initial condition, /9^(0) = Ig. Similarly, the stochastic complex time-dependent Hamiltonian is 
given by 



(26) 


n 


where the scalar noise components are again of zero mean and correlation: 


(^n(^)) = 0 , 

(^n(’2)^m(’2)) = dnmCni^Z Z )/^ . 


(27) 
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FIG. 1. The integration contours in the complex time plane z = t — it used in the various calculations. 
Integrating along z = t (green arrow) yields the absorption operator while the contour 2 ; = —it (red arrow) 
results in the equilibrium reduced matrix. The blue arrows characterize the independent contours needed 
to generate the emission operator at times s and s'. 

As with the absorption, the numerically exact emission operator is obtained after the stochastic 
average over the noise variables, E^(t)* = {fP{t — ih/3))^/{Z)^ where Z = Tt [p^ {—ihl3)]. 

It should be noted that the stochastic path integral expression in Eq. [25] represents a gener¬ 
alized form of the stochastic propagation. It contains two interesting limits that are represented 
schematically in Fig. [TJ For example, if the imaginary time component of the propagation is set 
to zero such that z = t, then one immediately recovers the absorption operator obtained above in 
Eq. [T7| Alternatively, if the real time component of the integration contour is set to zero, then one 
recovers the pure imaginary time evolution of the equilibrium reduced density matrix propagation 
explored in our earlier worki^ This leads to the interesting result that the emission operator at 
t = 0 is simply the exact equilibrium reduced density matrix. 


D. Computational considerations 

There are several points with regards to the stochastic formulation that should be emphasized. 
Firstly, a generalized stochastic approach to compute the real time dynamics of the entire reduced 
density matrix has recently been explored.— In that case, the presence of complex noise gen¬ 
erally leads to very slow convergence of the stochastic average as the length of the simulation 
increases. The approach presented here, and in particular Eq. [m represents a simplified version of 
those works, and thus directly inherits their numerical difficulties. However, the redeeming feature 
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of the present approach is that the decay time of the absorption and emission correlation functions 
is much shorter than the corresponding relaxation time of the pure real time dynamics. As a re¬ 
sult, brute force convergence of the stochastic path integral is generally possible with a reasonable 
number of Monte Carlo samples (10^ — 10®), although the low temperature regime can be more 
demanding. 

Secondly, there is a subtle difficulty that should be discussed with regards to the calculation of 
the emission spectrum. As seen from Fig. [H the propagations to times s — ihf5 and s' — ihf5 evolve 
along different contours in the complex time plane. The bath correlation functions evaluated along 
these two contours are different, and thus the emission operators E(s) and E(s^) require completely 
independent calculations. In principle, this should increase the cost of computing the emission 
spectrum by a factor of the number of time steps with respect to that of the absorption. However, 
this is not the case since the presence of the imaginary time component in the propagation greatly 
improves the convergence properties of the Monte Carlo calculation^^iS^ While the computational 
cost of the emission spectrum is more expensive, it is not prohibitive. 

Finally, the inclusion of static disorder in the absorption operator calculation is trivial, but less 
so for that of the emission operator. In the former case, the averages over the noise and disorder 
commute, and thus the two may be computed simultaneously. That is, the disorder-averaged 
absorption spectrum should incur practically no additional computational cost over that of the 
bare absorption spectrum. However, the presence of the partition function in the denominator of 
the emission operator demands that the average over the disorder must be computed independently 
of the average over the noise. As a result the disorder-averaged emission spectrum, although 
straightforward, may be quite costly from a computational perspective. 


III. NUMERICAL RESULTS 

Although the path integral formalism is valid for any spectral density, below we will focus on 
the standard Drude-Lorentz form so that benchmark results from the HEOM formalism can be 
obtained. In a forthcoming work, we will examine the influence of the spectral density on the 
spectra and energy transfer rates of the light harvesting system LH2.— The Drude spectrum is 
defined by 


Jiu) = 

LO^ + 7 


(28) 
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where 7 is the cutoff frequency and the reorganization energy A is defined such that, 

J{uj) 


A = 


1 / 

TT Jo 


duj 

0 ^ 


( 29 ) 


As is commonly assumed, we take the spectral densities for each of the independent baths to be 
equivalent and, unless otherwise specified, fix the reorganization energy at A = 200 cm“^ and the 
cutoff frequency to 7 = 53 cm“^ (10 ps“^). 


A. Emission spectra 

Before presenting the MCFT rates, we will first focus on the far-field spectra. The absorption 
and emission spectrum can be computed by combining the knowledge of the corresponding oper¬ 
ators defined above in Eqns. [ 8 ] and [9] with the respective transition dipole moment vectors of the 
chromophores (/x), by 




/ OO 

dt 6+*^^* ^ (ti ■ pt) (e* • Pn) lmn(0 
n^m 

poo 

E? (w) = / dt [si ■ ) (Ti • /I°) E°„(t) 


(30) 

(31) 


where e* is a unit vector characterizing the polarization of the incident radiation field that projects 
onto the dipole moment vector of each chromophore. As noted above, the path integral evaluation 
of the emission spectrum contains the absorption spectrum as a limiting case, and hence we will 
focus only on the former here. 

As a preliminary benchmark calculation to prove the efficacy of the path integral approach. 
Fig. [2] displays the stochastic path integral results for the emission spectrum in the unbiased two 
level system that was explored in papers I and II. The system Hamiltonian Hg = Vax with V = 200 
cm“^, leads to highly delocalized exciton states, and thus serves as an interesting test case to assess 
the validity of the MCFT formalism as well as that of approximate perturbative approaches. For 
simplicity, the dipole moment operators have been chosen to be equivalent for each site and in each 
direction such that firn = 1 in Eons. 1501 and IHTl Because of this choice, the spectra are determined 
from the simple sum of all of the elements of the emission operator (cf. Eg. [3T]) . In Fig. [2l the path 
integral results at 200 and 300 K are compared with the corresponding results from the standard 
HEOM approach shown with increasing number of Matsubara terms. The HEOM results are seen 
to eventually converge to the path integral results, although even for this relatively simple two level 
system, the hierarchy results are difficult to converge and require both a large number of hierarchy 
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FIG. 2. (a) The emission spectrum at T = 100, 200, and 300 K in a model two level system comparing 

the present stochastic path integral calculations with the corresponding HEOM results with 0 (red) and 2 
(blue) Matsubara terms. The number of hierarchy tiers required for convergence in each case is 12. The 
bath is defined by the reorganization energy A = 200 cm“^, and cutoff frequency, 7 = 53 cm“^. For clarity, 
the results at T = 200, and 300 K are vertically offset by 0.0075, and 0.015, respectively. At T = 100 K we 
cannot converge the standard HEOM and only the stochastic HEOM (green) formalism can produce reliable 
results, (b) The reorganization energy dependence of the emission spectrum at T = 100 K. The remaining 
parameters are the same as in (a). 


tiers as well as several Matsubara terms. At the lowest temperature of T = 100 K shown in Fig. [2l 
the standard hierarchy calculations can not be converged with respect to the number of Matsubara 
terms, and the sHEOM approach must be used to generate the numerically exact results.—*^ Note 
that at each temperature, the hierarchy results and present path integral results are in precise 
agreement. However, compared with the hierarchy calculations, the stochastic approach developed 
here is more straightforward both in terms of implementation and convergence. Additionally, since 
the stochastic formalism is a Monte Carlo method, it is trivially parallelized and free from the 
memory demands that plague other density matrix approaches such as the HEOM or QUAPI. In 
the case that further improvements to the computational efficiency of the emission path integral 
are necessary, a very useful and accurate approximation can be employed which is discussed in 
Appendix [B1 

As is readily seen, the spectra in Fig. [2^ are comprised of two peaks centered around the 
eigenstates of the total system Hamiltonian. While the intensity of the peak at positive frequencies 
is nearly independent of temperature, that of the low energy peak steadily decreases and vanishes 
at the lowest temperature shown of T = 100 K. This is in stark contrast to the behavior expected 
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from an isolated two level system where the emission spectrum can be computed analytically as 

2 

E{u:) = + 6ico - ei) , (32) 

i=l ^ ^ 

with the eigenstate energies, ei ^2 = + A^, V is the electronic coupling, A is the bias, and 

the eigenstate population. Pi = As seen from Eq. [32] the spectra are composed of two 

peaks centered at the eigenfrequencies of the system with intensities that are weighted by the 
respective Boltzmann populations of the two states. At low temperature the population localizes 
in the ground state, and the spectrum shifts to the red. In Fig. [J] the opposite occurs and a 
blue shift is clearly seen with decreasing temperature. This behavior is a result of the strong 
system bath coupling. To demonstrate this effect more clearly. Fig. [2)3 displays the reorganization 
energy dependence of the emission spectra at the lowest temperature shown in Fig. [2)i of T = 100 
K. Only at very weak coupling, does the spectra resemble that expected for the isolated system 
from Eq. [3^ with the emission dominated by the low energy eigenstate of the system. However, 
as the coupling increases in Fig. [2)), the weighting of the two peaks is redistributed towards the 
higher lying eigenstate resulting in a steady shift to the blue. As discussed in papers I and II, 
the equilibrium state of the total system and bath cannot be written in a factorized form as in 
Eq. [32] particularly when the temperature is low and the system-bath coupling large, as is the case 
here. This is the key feature that is responsible for the drastic failure of standard perturbative 
approximations to the emission spectra as well as the counterintuitive temperature dependence 
seen in Fig. [2)i. 

B. MCFT rates 

We next consider the multi-chromophoric energy transfer, where both the donor and acceptor 
complexes are comprised of the symmetric two level system analyzed in Fig. [2] That is, each 
complex is described by the system Hamiltonian, Hg = Vcr^, and the donor-acceptor couplings 
(Eq. [2]) are set to the constant value, = 10 cm“^. This weak coupling ensures that the 

perturbative MCFT formalism is valid and is also characteristic of many natural systems.— The 
energy transfer rates computed as a function of the system-bath coupling strength are displayed 
in Fig. [3)i. Although not shown, the rates from the HEOM formalism are in precise agreement 
with the present path integral results in the region where the former can be converged (up to 
A = 600 cm“^). As the transfer occurs between two symmetric systems, the transfer rates are 
monotonically decreasing functions of the system-bath coupling strength as would be expected 
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FIG. 3. (a) Energy transfer rates between two symmetric two level systems as a function of the reorgani¬ 

zation energy. The cutoff frequency, 7 = 53 cm“^ and the temperature is T = 300 K. The black dots, green 
lines, red triangles, and blue squares denote the results from the present stochastic path integral approach, 
the second order time-convolution (TC2) master equation^ the full cumulant expansion (FCE))^ii and the 
hybrid cumulant expansion (HCE))^ respectively, (b) The energy transfer rates as a function of the temper¬ 
ature. The reorganization energy is A = 200 cm“^ and the cutoff frequency, 7 = 53 cm“^. The results from 
the TC2 approach lie outside the scale of the graph for all temperatures shown. 

from a simple analysis based on the standard Forster theory. 

Also included in Fig. [3^, is a comparison of the exact energy transfer rates with many of the 
commonly used perturbative methods. The TC2 is the standard second order, time-convolution 
master equation previously exploredJ^ As it is based upon the approximation of weak system- 
bath coupling, its validity is rather limited, and is generally inapplicable to many interesting 
physical systems, such as light-harvesting complexes, where the system-environment couplings can 
not be considered as small. The full cumulant approximation (FCE) explored in paper I and 
Ref. 3 provides reliable results over a much larger region of the parameter space as compared 
with the TC2, although it too begins to break down at very large system-bath couplings and 
eventually produces unphysical negative rates. The failure of both the TC2 and the FCE lies in 
their inaccurate treatment of the correlated initial state. Clearly, a perturbative expansion around a 
factorized initial state is qualitatively incorrect at large-system bath coupling. In order to overcome 
this difficulty, paper II explored an expansion around the numerically exact equilibrium reduced 
density matrix, which can be straightforwardly obtained through imaginary time path integral 
techniques.— As seen in Fig. [3^ this HCE technique provides a uniformly reliable approximation to 
the energy transfer rate, even at very strong system bath couplings. 

Fig. [ 3)3 displays the temperature dependence of the MCFT rates. Qualitatively, the results 
follow the predictions of Marcus theory displaying a maximum as a function of temperature. How- 
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ever, the Marcus rate formula predicts a maximum at 2A//cb ~ 600 K which is considerably lower 
than that observed from the exact calculations. Additionally, Marcus theory predicts that the 
energy transfer rate should vanish as the temperature decreases to zero. This is clearly not borne 
out in the exact results as the MCFT rates decrease to a finite value at low temperature due to 
non-vanishing quantum fluctuations. As with the system-bath coupling dependence, the approxi¬ 
mate perturbative results are also included. It is seen that below room temperature, the accuracy 
of the FCE quickly degrades and eventually produces negative energy transfer rates. The results 
from the TC2 method are outside the scale of the graph at all temperatures, which is not entirely 
unexpected as the system-bath coupling strength here is of comparable magnitude to all the other 
system parameters. However, the hybrid method provides reliable results across the entire parame¬ 
ter range, and also captures the plateau in the rates at low temperature. In a forthcoming work, it 
is demonstrated that the HCE is capable of capturing the temperature and system-bath coupling 
dependence of the energy transfer rates between two LH2 complexes while the failure of the TC2 
and FCE is even more dramatic than seen here.— 


IV. CONCLUSIONS 

A stochastic path integral approach to compute the energy transfer rates between weakly cou¬ 
pled multi-chromophoric complexes has been presented. As a consequence of the MCET formalism, 
one also has immediate access to the exact steady-state far-held absorption and emission spectra 
of the respective donor and acceptor complexes. The calculations of the absorption and emission 
operators require only the straightforward numerical solution of a stochastic differential equation, 
and the only difficulty lies in the convergence of the Monte Carlo average. As opposed to many 
other numerically exact approaches, the method developed here is amenable to any form of the 
spectral density, and can readily treat the low temperature and strong coupling regimes. To our 
knowledge, the present path integral approach is the only method currently available that can 
accommodate such a broad range of system parameters in relatively large excitonic systems. 

The numerical results presented here provide a systematic analysis of the role of the temperature 
and system-bath coupling strength on the emission spectra and energy transfer rates in model multi- 
chromophoric systems. As seen in Fig. [3] the exact MCFT rates serve as a stringent benchmark 
for approximate analytic methods. Whereas the standard perturbative approaches often yield 
qualitatively incorrect results, the hybrid cumulant expansion (HCE) technique developed in paper 
II^ can provide uniformly reliable results for the energy transfer rates across a large range of the 
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physically accessible parameter space. 
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Appendix A: Noise sampling 

Sampling the complex noise required for the absorption and emission operators is not completely 
trivial. The main difficulty is that the bath correlation function must be reproduced by (^(t)^(tO) 
rather than the Hermitian form . To proceed, the correlation function can be split into 

its real and imaginary components, 

C{t)=Crit)+iCi{t) (Al) 

and the influence functional rewritten as 

F[cj] = exp- - f dt' f dt"V{a{t'))V{a{t")) [Cr(t'— t") + iCi{\t'— t"\)] . (A2) 

L Jo Jo J 

Hubbard-Stratonovich transformations are then applied to each term separately, leading to 
F[o-]=j D[C]u;^exp dt' dt”Cit')C~^{t' - ^ dt'V{o{t'))C,{t') 

X J V[i']w„ex.p y dt' J dt"i^{t')C~^{\t' -t"\)iy{t") + J dt'V{a{t'))v{t') , 

(A3) 

where and denote the respective normalization constants. Thus the noise characteristics are 

(C(t)) = 0 {vit)) = 0 

= Crit - t') {l'{t)l^{t')) = Ci{\t - t'\) , (A4) 

and the autocorrelation function of the combined process, ^{t) = ({t) + y/iv{t) is readily seen to 
reproduce the desired bath correlation function, C{t). 
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Numerically sampling the real noise governed by the correlation, Cr{t), is straightforward since 
this kernel is strictly positive semi-definite. Sampling such noise has been discussed in detail in 
our previous works(^>^ One simply filters white noise with a kernel computed from the Cholesky 
decomposition of the Toeplitz matrix constructed from Cr{t). 

Sampling the noise for the imaginary part of the correlation function is less straightforward. The 
kernel, Ci{t) is not positive definite since C'j(O) = 0, so that the Cholesky decompositio n ap proach 
is not applicable. To cope with this, we have employed the approach suggested in Refs 
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and 
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First an eigen decomposition of the correlation matrix, Ci^nm = Ci{\tn — tm\), is performed, and 
the diagonal eigenvalue matrix is sorted into a nonnegative (A"*") matrix and the remainder (A~), 
which are both of the same dimension as Cj, such that 


C* = U [A+ + A-] U'^ 


(A5) 


The positive components are sampled in the usual fashion by filtering the appropriate kernel with 
white noise, while the negative components are sampled by taking the absolute value of the negative 
eigenvalues followed by a rotation with the complex unit. The desired noise sequence is then given 
by 


i7= 


(A+) 


1/2 


+ i 


A-|V2 


n 


(A6) 


where represents a realization of independent white noise terms. Using the properties of white 
noise, it is readily seen that the autocorrelation of iy{t) faithfully reproduces the desired imaginary 
part of the bath correlation function. 


Appendix B: Approximate Emission 

A very accurate approximation to the emission operator can be made by simply ignoring the 
imaginary part of the bath correlation function in Eq. [22j This simplification generally reduces the 
number of Monte Carlo samples required to converge the stochastic path integral by at least an 
order of magnitude. For the purely real-time dynamics of the absorption operator, ignoring Ci{t) 
leads to an extended Haken-Strobl model which rarely provides satisfactory results. However, for 
emission operator, the real-time and imaginary-time dynamics are intertwined so that the analysis 
is more subtle. In this case, there are still non-unitary contributions to the dynamics even if the 
Hamiltonian is purely real due to the complex-time evolution. To better understand this seemingly 
drastic approximation, it is useful to analyze the complex-time bath correlation function. It is 
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FIG. 4. The exact two-level system emission spectra reproduced from Fig.[5](solid black) compared with the 
approximate emission spectra (dashed red) computed by ignoring the imaginary part of the bath correlation 
function. The parameters are identical to those in Fig. [5J 

readily seen that the bath correlation function evaluated along the imaginary time axis to z = —ih(3 
is a purely real quantity for any spectral density. This case corresponds to the equilibrium reduced 
density matrix which is a purely real quantity if the Hamiltonian is real. Thus for small real times 
-during which time the emission operator has often substantially decayed- the imaginary part 
of the correlation function is also negligibly small. This approximation is particularly accurate 
in the high temperature limit where the increasingly broad emission spectra are a result of the 
increasingly rapid decay of the emission operator. In summary, in short-time limit ignoring Ci{z) 
is a reasonable approximation. In Fig. Hlthe exact results for the emission spectrum of the two level 
system are reproduced from Fig. [2] along with the corresponding results from the approximation 
scheme discussed here where the imaginary part of the complex-time bath correlation function has 
been set to zero. Only at the lowest temperature of T = 100 K are there any significant differences 
between the exact and approximate emission spectra. In fact, comparison with Fig. [2] indicates 
that even at T = 300 K, the approximate emission spectrum is more accurate than the HEOM 
results computed without including Matsubara terms. 
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